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Abstract 

We apply Runge-Kutta methods to linear partial differential-algebraic equations of 
the form A ut(t, x) + B(u xx (t, x) + ru x (t, x)) + Cu(t, x) = f(t, x), where A, B,C € 
M. n,n and the matrix A is singular. We prove that under certain conditions the tem- 
poral convergence order of the fully discrete scheme depends on the time index of 
the partial differential-algebraic equation. In particular, fractional orders of conver- 
gence in time are encountered. Furthermore we show that the fully discrete scheme 
suffers an order reduction caused by the boundary conditions. Numerical examples 
confirm the theoretical results. 

Key words: Partial differential-algebraic equations, Coupled systems, Implicit 
Runge-Kutta methods, Convergence estimates 



1 Introduction 



In this paper we consider linear partial differential-algebraic equations (PDAEs) 
of the form 



A u t (t, x) + B (uxx(t, x) + ru x (t, x)) + C u(t, x) = f(t, x), (1) 
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where t G (to,t e ), x G Q = (—1,1) C M., A, B,C £ M n,n are constant matrices, 
rfl, u,/: [to, t e ] x f2 — >■ IR n . We are interested in cases where the matrix A 
is singular. The singularity of A leads to the differential- algebraic aspect. 
It will always be tacitly assumed that the exact solution is as often differen- 
tiable as the numerical analysis requires. 

In contrast to parabolic initial boundary value problems with regular matrices 
A and B, here we cannot prescribe initial and boundary values for all compo- 
nents of the solution vector, they have to fulfill certain consistency conditions. 
We consider one example: 

Example 1 Superconducting coil (see Marszalek/Trzaska, Campbell/ Marszalek 
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with x G (0,/), t > 0. ui(t,x) denotes the voltage, u 2 (t,x) denotes the diver- 
gence of the electric field strength within the coil. I is the length of the whole 
winding. L, C and D are further coil parameters. Transformation to a partial 
differential-algebraic system of first order in t yields 
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As initial conditions we choose 



wi(0,x) = I — 



CDEl\ CDE , , n . 

x -\ — — x , u 3 {0,x) 



6 



6/ 



and as boundary conditions 
ui(t,0) = u 2 (t,0) 



0, m(t,l) = E, u 2 (t,l) = CDE, 



where E is the energizing source voltage at the input of the coil. 
As the boundary values of ui and m 2 are constant, we get from the third 
and fourth equation of ID) that u 3 and u 4 fulfill homogeneous boundary con- 
ditions. From the initial condition of ui and the first equation we derive 
^2(0, x) = CD ] E x. With ^3(0, x) = and the third equation it follows wit(0, x) = 
and therefore Ui xxt (0, x) = 0. With the first equation this implies M2t(0, x) = 0, 
and with the fourth equation we get finally 1*4(0, x) = 0. 

Here we have chosen the prescribed initial and boundary values such that all 
initial and boundary values are compatible. 
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For further examples considering the determination of the initial and boundary 
values which cannot be prescribed see Lucht/S./Eichler-Liebenow [7]. 

In the following we assume that for the numerical computation all initial values 

u(to,x) = <p(x), a;GO, 

and all boundary values entering into the space discretization are known, 

Bu(t,x) = ip(t,x), x G dQ, t G [ioi^e], 

where we restrict ourselves to Dirichlet boundary conditions to simplify the 
presentation. 

Investigations of the convergence of Runge-Kutta methods applied to ab- 
stract parabolic differential equations can be found for example in Bren- 
ner /Crouzeix/Thomee PQ, Lubich/Ostermann [6] and Ostermann/Thalhammer 
[10] . The approach used there cannot be carried forward directly to the class 
of problems considered here because the matrix A is singular. 

This paper is organized as follows: In Section [2] we derive a semi-discrete 
system based on finite differences. The result is a method-of-lines-DAE (MOL- 
DAE). 

Section [3] is devoted to the Runge-Kutta approximation of the MOL-DAE. 
Under a regular transformation, the MOL-DAE of dimension nN is decoupled 
into N systems of dimension n, where iV denotes the number of grid points 
on the x-axis. Furthermore, a Weierstrass-Kronecker transformation is used to 
decouple each of these systems into an ODE-system and an algebraic system. 
We introduce the differential time index of the linear PDAE and give the 
Runge-Kutta approximation to these subsystems. 

In Section H] we prove the convergence of L-stable Runge-Kutta discretizations 
with constant step sizes. The attained order of convergence in time depends 
on the differential time index of the PDAE and on the boundary conditions 
(homogeneous or inhomogeneous) which enter into the space discretization. 
Numerical experiments are finally presented in Section |5j We illustrate our 
convergence results for the backward Euler method and the 3-stage Radau 
IIA method. 



2 Space discretization 

The discretization in space of problem (JT|) by means of finite-differences results 
in a differential-algebraic equation (MOL-DAE) 

MU = DU(t) + F(t), t <t< t e , (3) 
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where U(t) is an iVn-dimensional real vector consisting of approximations to 
u at the grid points. Here N denotes the number of grid points on the x-axis. 
The matrix M is given by M = 1^ ® A and the matrix D originates from 
the discretization of the differential operator B-^-% by second order difference- 
approximations, from the discretization of the differential operator rB-^ by 
second (5 = |) or first order difference- approximations (5 G [0, 1] \ {|}) and 
from the matrix C, i.e., D is given by 

D = -—P®B-I N ®C, 

where In is the iV- dimensional identity matrix, 

( -{2 - hr{\ - 28)) l + hr5 
, l + hr(S-l) -(2 - hr(l - 2S)) \ -\- IirS 



\ ... ... l + hr(5- 1) -(2-hr(l-5))J 



and h = j^L. denotes the constant grid size. The Nn- dimensional real vector 
F(t) arises from the right hand side / of ([T]) and the boundary values which 
enter into the discretization. 

We denote by Uh(t) the restriction of u(t,x) to the spatial grid and by oih{t) 
the space truncation error defined by 

a h (t) := MU h (t) - DU h (t) - F(t). (4) 

By Taylor expansion of the exact solution we get 

a h (t) = h p *(I N ®B)j h (t) with KWHoo < K, (5) 

where p x G {1, 2} is the order of approximation of the space discretization and 
K is a positive constant, i.e., 

a h {t) = 0(h Px ) as h->0. 

Furthermore, we can show that there exists a regular matrix Q with 

Q" 1 ^Q = diag{A 1 ,...,A JV }, (6) 

where 



^ 2-hr(l-2S) +c l + hr5 1 + hr(5 - 1) ^ jir 



h 2 h 2 V l + hr5 N + 1' 

In the discrete L2-norm we have 

maxiJIQIUIQ- 1 !!}^ (7) 
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with a positive constant C\ independent of h. Therefore, in the following this 
norm is used. 



3 Runge-Kutta approximations 



In order to numerically advance in time the solution of the MOL-DAE ([3]), we 
employ an s-stage Runge-Kutta method 



U$ +1 = U m +T y £a ij K^. 1 , MK% +1 = DU$ +1 + F(t m + qr), 1 < i < s, 



U m+ i = U m + T^2biK® +1 , 
i=i 

where Oy, bi, Ci G R are the coefficients of the method and r = ^=*a ^j me 
step size. 

For the investigation of the convergence of the method, it is useful to intro- 
duce the Runge-Kutta matrix 21 = (diy)lj=i an d the vector notation t s = 
(1,...,1) T GM S , b=(b u --- ,b s ) T . 

Then, with the Kronecker product, we obtain the compact scheme 



U m +i = U m + T\b ®I Nn )K m+ i, (8a) 
S m+1 = t s ®U m + T{%®I Nn )K m+u (8b) 
(I s ® M) K m+l = (I, ® D) S m+l + F(t m+1 ), m = 0, . . . , M e - 1, (8c) 

where Sm+i = f ) • • • j ^m+1 1 ) = j • • • ) 1 an d 

F(Wl) = ( F (tm + ClT) T , . . . , F (tm + C S T) 

By the regular transformation ([6]), the MOL-DAE ([3]) can be decoupled into 
N DAEs 

AU Qk {t) = D k U Qk (t) + F Qk (t), k = 1, . . . , N, (9) 
with D k = —X k B — C and 

(U Q1 (t) T , U QN (t) T ) T = (Q- 1 <g> In) U(t), 

(F Q1 (t) T , F QN (t) T ) T = {Q- 1 ® I n ) F(t). 

In the following we assume that the matrix pencil {D + XM}, X € C, is regular, 
which is equivalent to the regularity of all the matrix pencils {D k + XA}. 

Definition 2 Suppose that all matrix pencils {D k + XA}, k = 1, . . . , N, are 
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regular and have the same index v^t- Then the differential time index of the 
linear PDAE (QJ) is defined to be v^. 

According to Weierstrass and Kronecker there exist regular matrices P k and 
Qk with 



P k AQ k = diag{/ nfel , . . . , I nksk , N mkl N mklk }, 

PkDkQk = diag{i?fci, . . . , Rks k , Im kl , ■ ■ ■ , ^m fc!fc }> 

where 



(10a) 
(10b) 



Rki 



^ki 1 







Xki 1 







(Z (T n ki, n ki AT 

e ^ , J-y mki 



^ki J 



1 







1 




G C 



(10c) 



(see Hairer/ Wanner [5]), and for the differential time index of the PDAE it 
follows 

v dt = max{m fei : % = 1, . . . , l k }. 



Therefore, DAE (0 is decoupled into systems of the form 

U\ki(f) = RkiUikiit) + F 1M (t), I = 1, . . . , Sfc, 
N mM U 2 m (t) = U 2 ki(t)+F 2k i{t), 1 = 1,..., l k 



(11a) 
(lib) 



with 



and 



(u lk i(t) T , U lkSk (t) T , U 2k i(t) T , U 2k i k (t) 



Qk'u* 



Qk 



[Fiki(t) T , . . . , Fi kSk (t) T , F 2 ki(t) J , 
Similarly, DAE (jl]) can be transformed to 



F 2M k (t) 



P k F, 



k r Qk- 



Uhiki{t) — RkiUhikiit) + Fi k i(t) + a h i k i(t), I — 1, . . . , s k , 
N mkl U h2k i{t) = U h2k i{t) + F 2k i(t) + a h2k i(t), I = 1, . . . , l k . 



(12a) 
(12b) 



Runge-Kutta methods are invariant under the transformations and ffTUj) . 
Therefore, to analyze convergence it is sufficient to apply them to systems of 
the form (11 lap and (lllbl) . Application to filial) yields 



U\kl,m+l — UiM,m + t \ b (g) I nkl j Ki k ^ m+ i, (13a) 
S\kl,m+1 = %® U lk l,m + T (21 <g> I nkl ) Kikl, m+ 1, (13b) 

Kiki,m+i = (Is ® Rki) S lk i, m+ i + F lk i(t m+1 ), m = 0, . . . , M e - 1, (13c) 
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and to (lllbj) 



u. 



2kl,m+l 



u. 



2kl.m 



+ T[b l (8)7, 



rikl 



K 



2kl ,m+l ; 



S2M ,m+l h® U 2k i, m + T{Ql®I nkl ) K 2kl 

,m+l 5 



[I s <8> ^m w ) K2kl,m+1 — S2kl,m+1 + i*2jw(£m+l) 



0, . . • , M e 



(14a) 
(14b) 
(14c) 



Now we start our convergence investigations. 



4 Convergence estimates 

At first we introduce the global (space-time discretization) error e m+ i and the 
residual (space-time discretization) errors 5 m+ i and A m+1 at the time level 

t tm+l- 

Definition 3 The global error e m+ i at t m+ i is defined by 

e m+l := Uh{t m+ i) — U m+ i 

and the residual errors S m+ i, A m+ i are given by 

<Wi ■= U h (t m + r) - U h (t m ) - t (b T <g> I Nn ) K m +u (15a) 



A m+ i := S m+1 -l s ® U h (t m ) - t (21 ® I Nn ) K m+l , (15b) 
where S m+ i and K m +i are defined by the exact solution Uh(t) of the PDAE, 
i.e., 



S m+ i := (U h (t m + Cit) 1 , . . . U h (t m + c s t ' 1 



K m+ i := (Uh(t m + Cir) T , . . . U h (t m + c s t) 



T 



Definition 4 The discretization scheme |2]) is convergent of order (p x ,P*), if 
the global error satisfies 

||e m+ i|| = 0{h Px ) + 0(t p *) for (m + l)r = const., r, h -»• 0, 

whenever u(t, x) is sufficiently often differentiate. 

With the components eQfc )m+ i defined by 

( e Ql,m+l) ■ ■ ■ 1 e QN,m+l) := (Q 1 ® 7i) e m+l 

and (JTj) we obtain the estimate 



1 



N 



7T\ ^ ll e Qfc,m+i|| 2 < || e m+l|| < ClA n 2^ 

°1 \ fc=l \ k=l 



N 



h J2 ll e Qfc,m+l|| 2 - (16) 
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Letting a hlMjm+1 = [a h iki(t m + cir) T , . . . , a hlk i(t m + c s r) T j we get with 
Pg) and (dScJ) 



A' 



,m+l ,m+l 



(/ s (g) ifo) (S 

lkl,m+l ,m+l + OChlkl ,m+l • 



Using (|13bp . the transformed components e\k\, m = Uhiki{t m ) — °f the 

global discretization error and the transformed components A lfc ^ m+1 of (115bl) 
we obtain 

Sikl,m+l~ S\kl ,m+l -U-s &lkl,m lkl,m+l ~ K\kl,m+1 ) +Aifci im+ i. 

Combining the last two equations leads to 
G(rR kl ) (k 

lkl,m+l — K\kl,m+\) — (^s ® -Rfci) &lkl .in ,m+l+ Q; Wfe/,m+l) 

(17) 

where G(z) — I s — z%. 

Remark 5 For a (matrix-valued) junction f(z) : C — >• C m,n which is analytic 
in a neighbourhood of Ku, the matrix function f(Rki) with R k i given in MUc\) 
is defined by (see Golub/van Loan J^j) 



f(Rki) 



f (is \ J ni2 ( * ki> 

JixilK-^ki) i\ 



fili2 {^ki 



(n w -l)! 



' «i=l,...,m, i2=l,...,r 



In the following we assume that the Runge-Kutta method is A-stable and 
^(>ffei) < or < C 2 for all h G (0, ho] with a positive constant C 2 . Then 
for sufficiently small r the matrix G{rRki) is regular, and the Runge-Kutta 
system (Tl3l) has a unique solution. Using (I13ap . (TTTh and the transformed 
components Siki, m +i of ( jl5a[) we obtain the recursion 

Glkl,m+1 = R(TRkl) e lkl,m + L(t R k l) A lfe ; >m+1 + TJ(TRkl)(Xhlkl,m+l + ^ lfcZ,m+l 

(18) 

for the discretization error eifc; )m+ i, where we have used the abbreviations 
J( z ) = b T G(z)-\ R(z) = 1 + J(z)l s z, L(z) = J(z)z 



(R(z) equals the classical stability function of the Runge-Kutta method). 



Published in Applied Numerical Mathematics 53 (2005) no. 2-4, pp. 213-229, 



doi: 10.101 6/j.apnum . 2004. 08. 023 



Solving the recursion (TT8T) with eo = leads to 

m 

e lkl,m+l = lkl,m+X-i 
i=0 

m m 

+ r^2 R{rRki) l J (rRki)a h iki,m+i-i + E R{TRki) % hki,m+i-i- (19) 

j=0 i=0 

Now we assume that the Runge-Kutta method under consideration has (clas- 
sical) order p and stage order q (p > q). Then the simplifying conditions (see 
Hairer/ Wanner [5]) 

B(p): jZh4- l = \, k = l,...,p, 

i=i K 

s 1 
C(q) ■■ ^J )' 1 = T c t i = 1, . . • , s, k = l,...,q, 

j=l K 

are fulfilled. 

With a Taylor expansion of Uh(t m +Cjr) and Uh(t m +Cjr), j = 1, . . . , s, around 
t m up to the order p we obtain for the j-th component of the residual error 
A m+ i the equation 

P r 



r=q+l 



with c l = (c\, . . . , c\) T . Therefore, with r A iw,m+i = \r\ llkh 
and 

_ L(g) g - rgg- 1 ] 
rl J ~ 1 - R(z) 
the error equation (fT9j) can be written as 



m +l, • • • j r A s lfcZ,m+l 



e lkl,m+l — T YR(rRkiYJ(rRki) OChlkl,m+l-i ,m+l—i 
i=0 i=0 
m 

+ Y R ( rR kiy L ( rR kl)rAlkl,ra+l 

8=0 

m p r 

+ 2 R(rRkiY (/„„ - i2(ri?w)) £ -^(^)^iL(^-) • 

j=0 r=ij+l r ' 

V v ' 

= K 

(20) 



Remark 6 The junction W r (z) was introduced by Ostermann/ Roche JM/ to in- 
vestigate the convergence of Runge-Kutta methods for abstract scalar parabolic 
differential equations. 
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For the subsequent error estimate, the term k is transformed in the following 
manner: By exchanging the order of summation we get 

P f r f m_1 

r=<j+l '• \ i=0 

m 

-Y^RtrRMT-^W^TR^U^U) - R(rR kl r +1 W r (rR kl )U^ kl (t ) 



From this we obtain 



V T r I m-1 



* = E 7 E ^(r^O'-Wr^) (tf£?„ft) - u$ kl (t i+1 

r=q+l r - \ i=0 

+ W r (r^)^iL(*m) - i2(ri? w ) m+1 W r (T^)^iL(*« 
Therefore it holds 

P r / m-1 

K = E - - E R(TRki) m - l W r (TR kl ) / U^P(s)ds 

r=q+l T - \ i=0 I 

+ W r (TR kl )U$ kl (t m ) - R(TR kl r +1 W r (rR kl )U$ kl (t c 

A similar transformation can be found in Brenner/Crouzeix/Thomee [I]. 
Inserting this into f )20|) results in 



eifc/,m+l = 7" E -^( r RklYJ (t Rkl)othlkl,m+l-i + E ^)^lfc/,m+l-i 
i=0 i=0 
m 

+ 53 R{TRki) l L{TR kl )r Alkl ^ m+1 

i=0 

p r / m-1 

+ E — - E ^(^) m -w,(^) / u§${8)d* 

r=q+l r - \ i=0 (. 

+ W r (TR kl )U$ kl (t m ) - R(rR kl r +1 W r (rR kl )ui\li(to) 



Assuming that the Runge-Kutta matrix 21 is regular we can derive an analo- 
gous equation for the components e2ki,m+i of the transformed global discretiza- 
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tion error 



m m 

^2kl,m+l = T ^2 R(NklY J(Nki)Olh2kl,m,+l-i + E R (-^kl) % ^2kl,m+l-i 
i=0 i=0 
m 

+ ^ H(JV fc i)^(^M)^A2JH,m+l 

ti + l 



P r / m— 1 r 

E - - E ^(^) m -w.(iv w ) / t/fisfta) rf s 

r= g +l r - \ i=Q / 



+ W r {N kl )U^ kl {t m ) - R{N kl ) m+1 W r {N kl )U^ kl {h) 



with the abbreviations 



J(z) = b T (I s z - r2l) _1 , = 1 + tJ{z)%, 



1 - R{z) 

Finally, using we get for e Qk>m+1 = Q k {eJ khm+1 , ej ksk>m+1 , ej klm+1) eJ Mfcjm+1 ) T 
the equation 

s m 

e Qk , m +i =h Px T E <5fcdiag{. . . , R(rR kji y Jj(TR kjl ), . . . , 

j=l i=0 

m-i + c j T ) 

m 

+ ^Q fc diag{. . .,R(TR kjl )\ . . .,R(N mk]2 )\ . . . }Q k 1 5 Qk;m+l ^ 

i=0 
s m 

+ EE^ dia s(' • ^ R { TR kh) lL j{ TR kh), 

3=1 i=0 

R(N mk jL j (N mk . 2 ),...}Q k 1 r Aj Q,m+1 
+ E ^Qk\di &gk {...,W r (rR kjl ),...,W r (N mk .^ 

r=q+l T - \ 
m-1 ti+1 

- E / diag fc {. . . , R(rR kn r- l W r (rR kn ) } R(N rrikj2 ) m - l W r (N mk32 ), ...} 

2 = /. 

- diag fc {. . . , R{TR kh ) m+x W r {TR kjl ), R{N mk]2 ) m+1 W r {N mk]2 ), ...} 

QfuSlfa)). (21) 
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Now we can estimate the different terms in ( 121]) . For that purpose we assume 
in the following that the matrix norms 

||Q,diag{iV; fei , o, . . . , ojQ^l. . . , ||Q fe diag{o, . . . , N^JQ^W (22a) 

and 

\\Q k diag{N* nja ,o, o}P k B\\, . . . , ||Q fc diag{o, . . . , o, }P k B\\ (22b) 

are bounded for i = 0, . . . , maxjV^, : ji = 1, . . . , s k } — 1 and all h G 
(0, ho] , where o denotes a zero matrix. 

Because of the A-stability of the Runge-Kutta method and ^(xjyj < or 
l^fcjj < C*2 for all h G (0, fto] we have that WR^rRkj^W, \\ Jji^Rkj^W and 
||Lj(T.Rfc.n)|| are bounded for sufficiently small r. We assume further that 
Riit) ^ 1 for t G R \ {0} and lim ^ 1. Then ^(r-RjyJ exists and 

z— >— oo 

is bounded. Moreover, as it is shown in Ostermann/Roche [9], one has 

\\T r W r (TR kjl )\\ = 0(^{P.ff+2+«})||i^{°." nin W.9+2 + a-r}}|| 

with a 6 R, a > -1. 

Assuming that \xkn\ < £3(1 + |Afe|) one can show (cf. D.[3], the proof relies 
on the Mean Value Theorem and Abel's partial summation formula) that 

ii^rn \\u { Qi(t m )w=k^o(h^). 

Altogether the terms in (T2~T!) that originate from e\ki, m +i are of order 

0{h p *) + o( T min{p ' q+2+a} )h-h 1+2a . (23) 
With the Taylor expansion 

l hQk (tm-i + C 3 T) = g g k\{l-k)\ ^dt llhQk{tm) + °^ VAt ^^ 

the term 

s m 

a =h p * J2 T E Qkdiag{0, . . . , 0, R{N mk .J Jj{N mkh ), 0, . . . , }P fcJ B TftQfe (t m _ i + c,-r) 

j=l i=0 
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of eQk, m +i can be written as 

j=0 J- 1=0 ' k=0 li -\ b 

m , ... o 

J( z )g*) w (0), 0, ... , 0}P fc 5^ 7 , Qfc (t m ) 
i=0 °^ 

+ E 7rE^diag{0,...,0,M (^ +1 fl(z) < J(z)) W (0),0,...,0}fl k SO(7 J *- 1 ; 

j=0 J- i=0 

where (. . . denotes the j-th derivative w.r.t. z. 

For L-stable Runge-Kutta methods with regular coefficient matrix 21 we have 
R(oo) = 1 — 6 T 2l _1 ]l = and therefore R(0) = 0. If the matrix norms in (j 2 2 aft 
are bounded, then ||a|| is bounded if 

E TJ7T mt EH) 1 -" (r J+1 i^)V>)£ fc f } (0) = (24) 



£5 *!(*-*)! 



i=0 



for Z = 0, . . . , j - 1, j = 1, . . . , u dt - 1. 

The remaining terms in the equation (pTI) yield the classical order of 
the Runge-Kutta method applied to a linear DAE of index v^t with constant 
coefficients. Thus, altogether we have 

e Qk ,m+i = 0{h p *) + 0(r p ^) + o(r mm{p ' q+2+a} )h~h l+2a . 

N 

From (fl6|) it follows that we have to choose a such that J2 k 2 ( 1+2a > is bounded 

k=l 

for N — > oo. This implies a = — | — e, e > 0, and we have 

||e m+1 || = 0(h P *) + 0( r min{ P „ dt , g +1.25-e}^ 

If the derivatives of order (q + 1) w.r.t. the time of the boundary conditions 
that enter into the space discretization are homogeneous, i.e. 

d q+1 u 

B°—^- = 0, (25) 

(DEi) yields 

RklU$ kl (t) = Ufc£\t) - F lk iit) + a hlM (t), 
and instead of (1231) we obtain the order 



0{h Px ) + o{T min{p ' q+2+a} )h-h 2a -\ 
which implies a = \ — e, and therefore 

||e m+ i|| = <D{h p *) + o(r miD{p "^ q+2 - 25 - £} ). 
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Summarized, we have the following convergence result for smooth enough 
solutions u(t,x) of the PDAE (max-jjo + l,p + u dt — 1} times differentiable 
with respect to t in [t , t e ] and p x + 2 times differentiable with respect to x in 
[-1,1]): 

Theorem 7 Let the following assumptions be fulfilled for h — > (N — > oo) 
and k = 1, . . . , N : 

a) for the matrix pencils D k + XA there exist Weierstrass-Kronecker decom- 
positions according to /[TO]) , and the matrix norms in &2}) are bounded, 

b) 3?(x fcjl ) < or \xkji\ < C 2 for all h e (0, ho], 

c) l^kjA < C 3 (l + |A fe |) ; 

d) if v<it > 2 then (24\ ) is fulfilled for I = 0, . . . , j — 1, j — 1, . . . , v^t — 1. 

Furthermore let the Runge-Kutta method be of consistency order p, stage order 
q and L-stable (ifv^t = or v^t = 1 , it suffices A-stability with lim Riz) < 1) 

with a regular matrix 21 and R(it) ^ 1 for i Gl\ {0}. Let p Vdt be the classical 
order of the Runge-Kutta method applied to a linear DAE of index Vdt with 
constant coefficients. 

Then the discretization method |5J) converges for linear PDAEs after v dt time 
steps with the order (p x ,P*) in the discrete L2-norm in space and in the max- 
imum norm in time with 

!mm{p u , q + 1.25 — e}: inhomog. boundary conditions according to l[2b)) 
mm{pv dt , q + 2.25 — e}: homog. boundary conditions according to [2h^l 

and e > arbitrary small. 

Remark 8 (1) Stage order q > Vdt — 2 implies condition d) in Theorem 1 
for v dt = 3 or v dt = 4. 

(2) The assumptions on the Runge-Kutta method are fulfilled, e.g., for the 
Radau IIA and the Lobatto IIIC methods and in the case of v dt < 1 also 
for the implicit midpoint rule. 

(3) If p > p* + 1, then for L-stable Runge-Kutta methods with zL(l,z) 
bounded for 3t(z) < 0, the condition that Ii22a\) is bounded can be replaced 
by the boundedness of the matrix norms 

|| — Q k di&g{N l n , o, . . . , o}Q k 1 1| , . . . , || ^Q fc diag{o, . . . , o, N l n , o, . . . , o}Q k 

(26a) 

|| Q fe diag{o, . . . , o, Nl ki , o, . . . , o}Q- k 1 \\ , || Q fc diag{o, . . . , }Q k 1 1| . 

(26b) 
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(4) If we choose e = 0, then we get for p* < p- 



"at 



\e m+1 \\=0(h p *) + 0(J\lnh\TV*). 



Remark 9 For a given Runge-Kutta method, Theorem [7| can be specialized. 
E.g., if we take the implicit Euler method, the resulting BTCS method is con- 
vergent of time order 1 for arbitrary time index, if only the conditions a) and 
b) of Theorem^ are fulfilled. For the Radau II A methods with s > 2 stages we 
get 



min{s + 1.25 — e, 2s — 1}: = 0, 1, inhomog. b.c.s according to ( f^3]j 
mm{s + 2.25 — e, 2s — 1}: u^t = 0, 1, homog. b.c.s according to §Wt 
s + 2 - v dt : v dt >2 



as temporal order of convergence, provided that the assumptions a)-d) of The- 
orem [7| are fulfilled. 



5 Numerical examples 



The numerical examples given below illustrate our convergence results. For the 
time integration we use the backward Euler method and the code RADAU5, 
which is a variable step size implementation of the 3-stage Radau IIA method, 
see Hairer/Wanner [5j. The Euler and Radau IIA methods are of great impor- 
tance in applications. 

Example 10 The backward Euler method is given by the parameters 

s = l,2t=(l),6 = c=(l),p = g=l. 
We consider the linear PDAE 



f Q 1 




(° 





-A 




/ 


-1 


-1 


-A 




1 


Ut + 





-1 


-1 









-1 





u 


v° °y 




.° 





°y 




I 












=A 






=B 










=c 







a coupled system of two parabolic equations and one algebraic equation, with 
x G [—0.5,0.5], t G [0, 1]. The right-hand side, initial and Dirichlet boundary 
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values are chosen such that 



u(t, x) 



is the exact solution. It holds 



1 x(x — 1) sin(i) ^ 

x(x — 1) cos(i) 
\yx(x - l)(e* + t 5 ) J 



Du = —XhB — C 



11 1 + A fc 
1 + A fc A fc 
1 



With 



\ 



A fc + 1 





-1 
1 




-Afc — Ajt — 1 
-A fe - 1 
1 



/ i 



and Qk 



Afc + l 







1 




1 



V 







A fc +1 A fc +1 

1 



we obtain the Weierstrass-Kronecker decomposition 

I 



PkAQk 



^0 1 0^ 

1 
y0 0y 



PkDkQh 



1 
1 



Therefore, the PDAE has differential time index 3, and the assumptions |HP 
are fulfilled. Remark^ yields that the BTCS method is convergent after three 
steps of time order 1 . This is confirmed by the numerical experiment, Table U\ 
shows the observed order of convergence in time at (x = 1, t e = 1). The nota- 



O.It- 1 


2 2 


2 3 


2 4 


2 5 


2 6 


2 7 


O.lh- 1 














2 2 


0.81 


0.91 


0.96 


0.98 


0.99 


0.99 


2 3 


0.81 


0.91 


0.96 


0.98 


0.99 


0.99 


2 4 


0.81 


0.91 


0.96 


0.98 


0.99 


0.99 



Table 1 

Numerically observed order of convergence in the discrete L2-norm. 

tion of the first element 0.81 denotes the observed order when refining the grid 
from (h = 0.1/2 2 ,r = 0.1/2) to (h = 0.1/2 2 ,r = 0.1/2 2 ), i.e., 0.81 = log 2 £ 7 
where £ denotes the ratio of the error with [h = 0.1/2 2 ,r = 0.1/2) to the 
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error with (h = 0.1/2 2 ,r = 0.1/2 2 ). Furthermore, we see that a simultaneous 
refinement of h and r yields no order reduction. 

Example 11 We consider the 3-stage Radau II A method with consistency 
order p = 5 and stage order q = 3, and the linear PDAE 



^0 2 0^ 



1 -1 
1 -1 

=A 



( 
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°) 




f° 


0^ 










^xx H" 





-1 
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Olj 




=B 








=c 





u = f(t,x) 



with x G [—1, 1] and t 6 [0, 1]. This example shows the dependence of the time 
order on the boundary values. 

1. We choose the right-hand side such that 

u(t,x) = (x 2 e _i , x 2 e~2 t , x 2 sin tj 
is the exact solution. Then we have inhomogeneous boundary values 

i(t, =f1) = (e~*, e~2*, sint) T . 



u 



Furthermore it holds 

^XkO ^ 
D k = 1 

A fe - 1 



with 



\ 



^k~Vk 
Xk+Vk 





P k AQ k 



J 






V 



i i 

1-A fc A fc -l/ 



) Qh 



^1 0^ 

1 
yO Oy 



4A fc 
(Xk+Vk)Vk 

Afc_ 

Vk 



PkDkQk 



2X k 









2A fc 











-Afe+r7 fe 
1 



4A fc 







{Xk~Vk)Vk 
Afe o 

Vk 

A 



Vk 



\ r, k (X k -l) V k(Xk-l) J 



(\h + 8| > \ for N > 3, i.e. h< \). 

The PDAE has therefore differential time index 1, and the conditions 
a)-c) of Theorem^ are fulfilled which yields convergence of time order 
4.25 — e. This is confirmed by the numerical experiment, see Tabled 
2. If instead the right-hand side is chosen such that 

u(t,x) = [{x 2 - l)e~* + t 3 cos 2 x,x 2 e~^\ [x 2 - 1) sint + t 3 sin 2 x) 
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1 T -1 
U.J-/ 


2 1 


? 2 

Z, 


o3 










2 3 


4.27 


4.28 


4.30 


2 4 


4.26 


4.26 


4.26 


2 5 


4.26 


4.26 


4.25 


2 6 


4.26 


4.26 


4.25 



Table 2 

Numerically observed order of convergence in the discrete Z/2-norm for inhomoge- 
neous boundary values. 

is the exact solution, then we have inhomogeneous boundary values 
u(t,Tl) = (t 3 cos 2 l,e"^,t 3 sin 2 l) T , 

but the derivatives ofBu(t,x) of order 4 w.r.t. the time vanish. Therefore, 
we obtain the convergence order 5 in time, see Tabled 
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5.00 


5.00 



Table 3 

Numerically observed order of convergence in the discrete L2-norm for inhomoge- 
neous boundary values where the derivatives of Bu of order 4 w.r.t. the time vanish. 



Example 12 We consider the 3-stage Radau IIA method and the linear PDAE 
([IP describing the superconducting coil. It holds 
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with 
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The coil PDAE has therefore differential time index 2, and the conditions of 
Theorem \7\ ( with the matrix norms Ii22a\) replaced by l[2b}) ) are fulfilled which 
yields an order of convergence in time of 3. 

This is confirmed by the numerical experiment, see Table^ 
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2 2 


3.00 


3.00 


3.00 


2 3 


3.00 


3.00 


3.00 


2 4 


3.00 


3.00 


3.00 



Table 4 

Numerically observed order of convergence in the discrete L2-norm for the coil 
PDAE. 



6 Conclusion 



The attention has here been restricted to a class of linear partial differential- 
algebraic equations. We have given convergence results in dependence on the 
type of boundary values and the time index. When the error is measured in 
the discrete L 2 -norm over the whole domain, the convergence order in time of 
the Runge-Kutta method for a smooth solution is in general non-integer and 
smaller than the order expected for differential-algebraic equations of the same 
index. Some numerical examples were presented and confirm the theoretical 
convergence results. 

The extension of the analysis to the case of space d dimensional linear partial 
differential-algebraic equations of the form 

d 

Au t (t,x) + J2 B i ( u xiXi{t,x) +nu Xi (t,x)) +Cu(t,x) = f{t,x) 

i=l 

with x = (xi, ■ • • ,Xd) T and a cuboid as domain is possible, see D. [3], but 
becomes rather technical and offers no new insight. Furthermore, the con- 
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sideration of periodic boundary values is also possible. Here we could show 
as temporal convergence order the order of an ordinary differential-algebraic 
equation. In the case of Neumann boundary conditions, the temporal conver- 
gence order lies in between the order obtained for Dirichlet- and the order 
obtained for periodic boundary conditions. 

Future work in this area will be concerned with convergence investigations for 
semi-linear partial differential-algebraic equations. 

Acknowledgements 

The authors are very grateful to the referee for his comments and fruitful 
suggestions. 



References 

[1] P. Brenner, M. Crouzeix, V. Thomee: Single step methods for inhomogeneous 
linear differential equations in Banach space. R.A.I.R.O. Anal. Numer. 16 (1982) 
5-26. 

[2] S. L. Campbell, W. Marszalek: The Index of an Infinite Dimensional Implicit 
System. Mathematical and Computer Modelling of Dynamical Systems 5 (1999) 
18-42. 

[3] K. Debrabant: Numerische Behandlung linearer and semilinearer partieller 
differentiell-algebraischer Systeme mit Runge-Kutta-Verfahren. Dissertation, 
Martm-Luther-Universitat Halle- Wittenberg, 2004. 

[4] G. H. Golub, C. F. van Loan: Matrix Computations. Third Edition. The John 
Hopkins University Press. Baltimore, London, 1996. 

[5] E. Hairer, G. Wanner: Solving Ordinary Differential Equations II. Stiff and 
Differential - Algebraic Problems. Springer- Verlag Berlin. Heidelberg, 1996. 

[6] Ch. Lubich, A.Ostermann: Runge-Kutta Methods for parabolic equations and 
convolution quadrature. Math. Comp. 60 (1993) 105-131. 

[7] W. Lucht, K. Strehmel and C. Eichler-Liebenow: Indexes and special 
discretization methods for linear partial differential algebraic equations. BIT 
39 (1999) 484-512. 

[8] W. Marszalek, Z. W. Trzaska: Analysis of implicit hyperbolic multivariable 
systems. Appl. Math. Modelling 19 (1995) 400-410. 

[9] A. Ostermann, M. Roche: Runge-Kutta Methods for Partial Differential 
Equations and Fractional Order of Convergence. Math. Comp. 59 (1992) 403- 
420. 

[10] A. Ostermann, M. Thalhammer: Convergence of Runge-Kutta methods for 
nonlinear parabolic equations, Appl. Numer. Math. 42 (2002) 367-380. 



Published in Applied Numerical Mathematics 53 (2005) no. 2-4, pp. 213-229, 



doi: 10.101 6/j.apnum.2004. 08. 023 



